package scu.maqiang.eigenvalue;

import java.util.Arrays;

import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3H81;
import scu.maqiang.mesh.Mesh3H8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.EigenSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class ThreeScaleSteklovEVFine3DH {

	public static void main(String[] args) {
		String dir = "SteklovEigenValue\\";
		Mesh3H8 mesh3D = new Mesh3H8().cube(64, 64, 64);
		//mesh3D.displayNodes();
		//mesh3D.displayElements();
		//mesh3D.displayBoundarys();
		int n = 8;
		double eps = 1.0/n;
		ScalarFunc chFun = (xy, label, param) -> {
			double xCell = xy[0] * n - (int)(xy[0] * n);
			double yCell = xy[1] * n - (int)(xy[1] * n);
			double zCell = xy[2] * n - (int)(xy[2] * n);
			return (Math.abs(xCell- 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25 && Math.abs(zCell - 0.5) < 0.25)? 1.0:-1.0;	
		};
		mesh3D.removeElements(chFun, null, 8);
//		mesh3D.toTecplot(dir + "PerforatedDomain3DH.dat");
		
		
		FES3H81 fs = new FES3H81(mesh3D);
		int dofHeat = fs.GetNdof();
		SRMatrix A = new SRMatrix(dofHeat);
		SRMatrix B = new SRMatrix(dofHeat);
		fs.assembleStiff(new double[] {1.0}, BVPType.COMMON, A);
		fs.assembleHeatConvectSurface(new double[] {1.0}, BVPType.CONSIST_MASS, B, 8);
		double[] x = new double[dofHeat];
		Arrays.fill(x, 1.0);
		fs.applyBCZero(A, x, Direct.X, BCType.RCE, 1, 2, 3, 4, 5, 6);
		EigenSSolver solver = new EigenSSolver();

		int N = 10;
		double[] D = new double[N];
		double[][] V = new double[N][dofHeat];
		solver.inversePowerTran_PCGSSOR(A, B, x, 2.0, N, D, V);
		System.out.println(MVO.toString(D));
		mesh3D.toTecplot(dir + "SteklovEV3DH64_8.dat", V);
		
	}

}
